Semiempirical Two-Dimensional Model of the Bipolar Resistive Switching Process in Si-NCs/SiO2 Multilayers

In this work, the SET and RESET processes of bipolar resistive switching memories with silicon nanocrystals (Si-NCs) embedded in an oxide matrix is simulated by a stochastic model. This model is based on the estimation of two-dimensional oxygen vacancy configurations and their relationship with the resistive state. The simulation data are compared with the experimental current-voltage data of Si-NCs/SiO2 multilayer-based memristor devices. Devices with 1 and 3 Si-NCs/SiO2 bilayers were analyzed. The Si-NCs are assumed as agglomerates of fixed oxygen vacancies, which promote the formation of conductive filaments (CFs) through the multilayer according to the simulations. In fact, an intermediate resistive state was observed in the forming process (experimental and simulated) of the 3-BL device, which is explained by the preferential generation of oxygen vacancies in the sites that form the complete CFs, through Si-NCs.


Introduction
In recent years, the study of resistive switching (RS) memories (RSM) has gained importance because of its application in non-volatile memories. Some features of RSM, compared to devices based on floating gate memories, are low electrical power consumption during writing, small size, fast operation [1] and the crossbar architecture of RSM arrays that enables less fabrication steps [2]. Nevertheless, multiple resistive states, memory window, repeatability and endurance are some other characteristics that need to be improved for RS devices. Therefore, these challenges need to be studied through the modeling of the physicochemical phenomena [3].
The typical structures for RSM are based on metal-oxide-metal (MOM) capacitive structure or metal-oxide-semiconductor (MOS)-like devices. The application of an electrical potential to these RSM produce a soft electrical breakdown that allows the formation of some conductive filaments (CF) in the RS materials, switching (FORMING, SET) the device from a high resistive state (HRS, OFF) to a low resistive state (LRS, ON). Once the ON state is activated, the RSM can return to the OFF state in the same polarity (unipolar RS) or by changing the polarity (bipolar RS) if the CF is broken [4].
Different materials such as HfO 2 , ZnO, TiO 2 , SiN x − Pt, SnO 2 , NbO x or Hf 0.5 Zr 0.5 O 2 have been reported as active layers in RSM devices [5][6][7][8][9][10][11], but the implementation of materials used in the actual semiconductor industry, i.e., Si or SiO 2 , needs to be investigated. In 2010, Yao et al. reported the RS properties of a SiO x -based device [12] and it was related to the deoxidation-oxidation process of oxygen vacancies that compound the CFs, which results in a change of the valence state of silicon atoms [13]. RS devices with the formation of a metallic filament by the diffusion of metallic ions from electrodes have been also reported in Cu/SiO x /W devices [14]. Devices with the deoxidation-oxidation process are known as valence change memory (VCM) while those with diffusion of metallic ions are known as electrochemical metallization memory (ECM) [15]. Some works have demonstrated RS properties in SiO x (x < 2) films, where the CF was formed by the connection of Sinanocrystals (Si-NCs) embedded in the SiO x film, corresponding to VCM [16,17].
In previous works, we reported the experimental RS properties of Si-NCs through the use of Si/SiO 2 multilayer (ML) structures where the Si-NCs act as conductive nodes that control the CF formation [18,19]. In a recent study on the bipolar RS behavior for Si-NCs/SiO 2 ML, the following was observed: if the FORMING process is obtained at reverse bias (positive voltage to the upper electrode), the RS shows an asymmetric currentvoltage (I-V) curve, while a symmetric I-V curve is obtained for a forward bias (negative voltage to the upper electrode), the latter being similar to MOM devices [20]. The RS behavior obtained in this kind of device is related to the composition and structure of the Si/SiO 2 ML. X-ray photoelectron spectroscopy (XPS) and high resolution transmission electron microscopy (HRTEM) revealed the formation of Si-NCs surrounded by SiO x (x < 2) layers. The analysis of the XPS-Si2p spectra measured in the SiO x films that separate the Si-NCs shows a pronounced presence of Si-suboxides, predominantly in the form of the O 3 ≡ Si − Si species (Si 3+ ), which has been reported as a sub-oxide that favors the field-induced formation of oxygen vacancies (V O ) due to its asymmetrical nature [18,19]. The presence of Si-NCs and the different Si n+ suboxides (n = 1-3), which are oxygen vacancies, influence the formation and annihilation of the CF and therefore the SET/RESET processes. Therefore, the simulation of the electrical properties of Si-NCs with these Sisuboxides, considered as oxygen vacancies, within the SiO x films is required to understand the different FORMING, SET and RESET processes observed in the RS behavior.
There are models that can explain the RS phenomena in different active materials. From the computationally complex to simple, the RS models can be classified as the following [4]: (i) Ab initio methods that model the ion transport during the formation of CFs [21] and the calculation of the energy states properties due to V O [22]; (ii) finite element methods that can use the properties obtained from Ab initio simulations. These models can simulate the generation and recombination between oxygen ions (O ions ) and V O by differential equations and their role for the electrical current [23][24][25]; (iii) Semiempirical models that simplify the computational complexity of the Ab initio methods by physical approximation of the formation of CFs and charge transport through it. These models use physical and chemical properties of the materials that could be obtained from method (i) and from experimental results of RS devices [6,[26][27][28]; and (iv) Compact models that can use the memristor theory [29] or conceptual simplifications to simulate and describe the RS in electrical circuits. These models take ideas like Leon Chua's memristor relations [30], cylindrical CFs [31,32], hourglass quantum-point contact [33] or random circuit breaker network [34]. The (i) and (ii) models can require specialized software and their accuracy depends on the computational power and the relatively large time of calculation. In the case of models (iii) and (iv), the reduction of computational complexity is an important advantage and they have demonstrated a good approximation to experimental data [4]. The compact models (model (iv)) simulate individual CFs or the relation between electrical charge and magnetic flux to implement them in circuit simulation; nevertheless, these approximations do not allow for the consideration of specific features of oxide materials such as nanoclusters, which can act as charge traps.

Materials and Methods
In this work, we use a semiempirical two-dimensional (2-D) model based on the VCM in relation with the deoxidation-oxidation hypothesis [26] to emulate the bipolar RS properties of Si-NCs/SiO 2 ML structures where Si-NCs are assumed to be fixed V O that cannot be recombined with O ions . These devices are formed with Si/SiO 2 bilayers deposited by RF-Sputtering on p-type silicon substrate (1-5 Ω-cm) with Al electrodes (Figure 1a). The presence of the Si-NCs/SiO 2 ML structure is confirmed by the HRTEM micrograph, shown in Figure 1b. Therefore, the simulation of these devices at different values of voltage sweeps allowed us to obtain, under this assumption, the position of V O generated in a 2-D mesh. These oxygen vacancy configurations (V O C) can be used to estimate the resistive state based on the length and number of CFs formed during bipolar RS [35].

Bipolar Resistive Switching Based on Deoxidation-Oxidation Model
From the rigid ion-point model for ionic crystals [36], and since the Si-O bonds can have an ionic nature [37], the drift velocity v of O ions at high electric fields in oxide thin films can be defined as follows: where a is the mesh size of the 2-D simulation, and it is assumed to be the distance between the neighboring barriers for the drift of O ions , 1/t 0 is the vibration frequency of O ions , E Om is the height of the potential barriers during the migration of O ions , k is the Boltzmann constant, T J is the temperature due to the Joule heating, q is the fundamental charge, ϕ dri f t is an enhancement coefficient for the drift of O ions , and F H is the electric field assumed as homogeneous in the active layer. The direction of F H depends on the electrical polarity and it is normal to the surface of the substrate with a magnitude equal to |V L − V R |/L where L is the thickness of the whole active layer. Therefore, the migration of O ions is opposite to the F H . The temperature T J is assumed to be T J = T r + |(V L − V R )I|R th , with T r as the room temperature, (V L − V R ) the potential difference between the electrodes, I the device current and R th the thermal resistance of the CFs [38]. Figure 2 shows the V O Cs due to the deoxidation-oxidation process to emulate the bipolar RS process of the Si-NCs/SiO 2 -based RSM. Due to the electric field F, the O ions (gray circle) are pulled up from the equilibrium position in the mesh, drifting at velocity v and piling them up at the substrate-oxide interface, as shown in Figure 2a. The sites left by these oxygen atoms generate V O (orange sites), decreasing the valence state of near neighboring silicon atoms as well as reducing the x value of SiO x and then increasing their metallic behavior [13]. These chains of V O can connect the electrodes through Si-NCs (blue sites) and contribute to obtain the LRS (Figure 2a). Since the O ions are accumulated near the substrate, when an opposite differential potential is applied, these ions return to the SiO x and some V O s on the left side of the active layer are recombined allowing that the device returns to the HRS (Figure 2b) [39]. For a time t, the generation probability of V O (P G ij ) and the recombination probability between O ions and V O (P R ij ) on each 2-D mesh site (i, j) are [26]: where E Oe is the migration barrier of the equilibrium position of oxygen ions in the lattice, F nH ij is the non-homogeneous electric field at each mesh site (i, j) ( Figure 2), γ is an enhancement coefficient with a different value for the SET and RESET process, β R is a constant coefficient for the recombination process, f i depends on the distance between the mesh site (i, j) and the substrate or left electrode, i.e., x i = a × i, and the product vt is the distance that an O ions travels and L O is the decaying length of the O ions concentration in the left side of the active layer. Since the linear model of T J is related to the temperature of CFs [40], T J affects the velocity of O ions (v), defined in (1), and the recombination probability, defined in (3). On the other hand, in this model, the generation probability defined by (2) is not affected by the Joule heating. Therefore, only the temperature T r is considered [23].
F nH ij and f i are defined as [26]: where V L is the voltage applied to the Si-substrate, which is assumed to be equal to the voltage of the bottom aluminum electrode, V R is the voltage of the top aluminum electrode (Figure 1a), δ Vo,ij is the Kronecker's delta that can determine if the site (i, j) is a V O .
As we can see in Figure 3a, the stochastic process consists in the following: at each potential difference (V L − V R ) and current I, if the site (i, j) is not a fixed V O or Si-NCs site (blue circles in Figure 2), the probability P G ij from (2), P R ij from (3) and a random number Rnd ij are calculated. If this site is a V O and the P R ij is greater than Rnd ij , the site passes to non-V O , which indicates it is not a V O . If the site is a non-V O and the P G ij is greater than Rnd ij , the site passes to a V O . In some cases, the RSM devices require a current compliance (I compliance ) to ensure reversible RS and to avoid severe damage in the active material [41]. The control of the current compliance is carried out at a new (V L − V R ) with temperature T J (Joule effect from previous values of current and voltage), as follows ( Figure 3b): a copy of V O C is saved in V O C copy ; then a new V O C, due to the initial time t = t init used when I < I compliance (Equations (2), (3) and (5)), is obtained by the stochastic process. Then, the current I is calculated for this new V O C. If the I is greater than I compliance , the V O C copy is recovered and t init is reduced by a factor of 1/1.1 y where y increases in one unit. The same algorithm is repeated until I ∼ = I compliance or that t init be reduced to a factor of 1/1.1 100 , assuming the V O C does not change and that the I = I compliance .

Oxygen Vacancy Configurations to Obtain the Resistive States
The resistive state, in bipolar mode, can be estimated using the Mott hopping model [42] for each different V O C obtained at the potential difference (V L − V R ) by calculating the hopping rate of the CF for the row j (R j ) ( Figure 2): where a 0 is the attenuation length of the electron wave function. The largest CF weighs the most since the V O chain will have the smallest distance to the left electrode ( Figure 2). Considering the M rows, the resistive state N S of each V O C is: The estimation of N S is included in the step "Calculate current I" in the flowchart of

Conduction Mechanisms at Different Resistive States
Different works have reported that the conduction mechanism depends on the resistive state, the RS material and the electrodes [43]. Based on the linear pre-factors for different resistive states proposed by Hu et al. [6], we use two factors, f HRS and f LRS that depend on N S and that multiply the current from different conduction mechanisms at HRS and LRS. These factors are defined by the piecewise functions shown in Figure 4b.
For the HRS, the charge transport can be due to the Poole-Frenkel (P-F) mechanism, which accounts for the reduction of potential wells at high electric field F H on the thin oxide layer [44], which is defined as [45]: where the constant β = q 3 /πε with ε as the dielectric permittivity, µ the electron mobility, N C the effective density of states in the conduction band, and φ t the potential of the trap levels below the conduction band. K HRS is an enhancement factor for the P-F mechanism. For the LRS, the space charge limited current with Frenkel effect (SCLC-F) was used, which takes into account the high F H on the thin oxide layer [46]: where K LRS is an enhancement factor for the SCLC-F. K HRS and K LRS are constant values that must be determined by experimental results because the 2-D model tries to estimate the behavior of the 3-D real devices, and N S is a dynamic value that is affected by the two dimensional V O C approximation. Finally, the current I at each voltage (V L − V R ) with the top electrode area A is given by:

Simulation and Results
The bipolar RS behavior of SiO 2 /Si-NCs bilayers (BLs) (6 nm thick/6 nm diameter) with an upper layer of SiO 2 (10 nm thick) was studied using MOS-like devices, as schematized in Figure 1a. Two RSMs were simulated, devices with one and three BLs (labeled as 1-BL and 3-BL, respectively). For our simulation, we use the values given in Table 1.   The high F nH ij defined in (4) promotes the formation of CFs near the Si-NCs because they are assumed to be fixed groups of V O increasing P G ij (defined in (2)). In other words, the presence of fixed V O , which acts as a conductive material, reduces the equivalent distance between electrodes inducing a high electric field. As observed in Figures 5b and 6b, a larger quantity of CFs are formed in the 1-BL device than the 3-BL device since a higher electric field is obtained for the thinner oxide of 1-BL device (Table 1). In these devices, the CFs at LRS can be preferentially formed in the rows where a greater number of fixed V O s exist; that is, where the Si-NCs are present as indicated by the arrows in Figures 5b and 6b. In fact, experimental evidence of this effect has been reported before [12].
The simulated N S , as well as the experimental and simulated I-V curves obtained for 1-BL-and 3-BL-based memristors are shown in Figure 7. A detailed explanation of FORMING, RESET and SET processes are presented below.

FORMING and SET Process
The V FORMI NG is modulated by the initial number of V O s (V O init) randomly placed within the active material. The 3-BL device is formed by two additional SiO 2 /Si-NCs bilayers compared to the 1-BL device. Therefore, a higher concentration of defects and V O need to be considered for this device [47]. The V Oinit value was 40 for 1-BL device whereas V Oinit = 400 was considered for the 3-BL to adjust the experimental V FORMI NG ( Table 1).
As we can see in the Figure 7b,d, the simulation of the FORMING process in both devices adjust to their experimental I-V data. The experimental FORMING process in the 1-BL device is observed as an abrupt change of current at V FORMI NG = 3.9 V while the simulated one occurs at 3.8 V (Figure 7a,b). For the 3-BL device, the FORMING process exhibits an intermediate resistive state (IRS) level, Figure 7d, until the current increases to the highest current value at V FORMI NG = 4.9 V (Figure 7c,d). As we can see, the simulated data fit well to the experimental data at the same V FORMI NG . The experimental IRS is obtained when the current jumps from 1.61 pA at 4.4 V to 627.07 pA at 4.5 V. Figure 8 shows the simulation of the V O C in this current jump. As we can see, new V O s (green dots) are generated as the simulated voltage increases from 4.3 V (indicated by A in Figure 7d) to 5.0 V (indicated by B in Figure 7d). Then, the IRS level can be explained as a result of new V O , preferentially generated at rows j = 25 and 34, as indicated by the arrows, which will form the complete CF to obtain the LRS (Figure 6b). The simulation of the 3-BL device (Figure 7c,d) also shows the change from the HRS state to LRS (FORMING). This process occurs when the voltage sweep returns to 0V (sweep 2) at V FORMI NG = 4.9 V, while the experimental FORMING is obtained for the voltage sweep from 0 V to 5.1 V (sweep 1) at V FORMI NG = 4.9 V (Figure 7d). This difference between the experimental and simulated FORMING process can be explained by the stochastic nature of the RS phenomena [48].
The simulation of the SET process for the 1-BL and 3-BL devices fits well to the experimental I-V data, as observed in Figure 7b,d, respectively. For the 1-BL device, the experimental and simulated V SET is about 2.7 V while for the 3-BL device, the simulated V SET is about 3.2 V, which is near the experimental V SET at 3.3 V.

RESET Process
Once the devices are in the ON state (LRS), an opposite voltage sweep produces the RESET (OFF) process. However, the simulations of this process for both 1-BL and 3-BL devices do not fit to the experimental RESET. Experimental V RESET values were −2.4 V and −2.7 V for the 1-BL and 3-BL devices, respectively. Meanwhile, simulated V RESET values were −3.1 V and −3.8 V for the 1-BL and 3-BL devices, respectively. This difference can be related to the stochastic process of the generation and recombination of V O during the RESET event. In fact, the RESET process has been reported to be dependent on thermal noise and nonlinear relaxation phenomena [49], which is not included in the present model.
As we can see in Figure 7b, the experimental RESET process of the 1-BL device was obtained when the voltage is swept from −4.0 V to 0 V (pink zone). The simulation of this process also shows the RESET at the same voltage sweep (green zone). Figure 9 shows the evolution of the V O Cs simulated for a voltage cycle from −3.2 V to −4.0 V (sweep 3) to −3.1 V (sweep 4) to explain the RESET process for the 1-BL device. When the voltage increases from −3.2 V to −4.0 V (Figure 9a (1) increases due to a T J = 443.13 K. Therefore, the distance that O ions can travel is about vt ∼ = 1446a, which is higher than the thickness of the active layer L = 45a (Figure 9b). The V O C remains unchanged as the voltage decreases from −4.0 V to −3.2 V (Figure 9b,c). However, the T J value decreases from 443.13 K to 336.4 K resulting in the decreasing of the distance that O ions travel to vt = 0.43a. Consequently, the O ions do not recombine with those V O s far away from the left electrode. As the voltage reduces from −3.2 V to −3.1 V (Figure 9c,d), an abrupt current drop is obtained from 658.17 nA to 3.44 pA (RESET). This strong current drop is related to an increased temperature T J from 336.4 K to 465.49 K, which results in a vt = 39a that allows the V O s to recombine, producing the RS from LRS a HRS (Figure 9d). As observed, the simulation of the RESET process depends on v of the O ions and the Joule heating since both control the recombination probability, as defined by (3). In fact, Figure 10 shows the exponential factor of the recombination probability for the 1-BL device.
A maximum value is observed at (V L − V R ) = −3.2 V as the voltage sweep increases from −2.5 V to −4.0 V. Nevertheless, the distance vt ∼ = 5a that O ions travel is not enough to obtain the HRS (Figure 9a). A second maximum value of the exponential factor is obtained at about −3.1 V when the voltage sweep reduces from −4.0 V to −2.5 V. This second maximum produces the RESET process due to the distance vt ∼ = 39a, which is similar to the thickness of the active material L = 45a. For (V L − V R ) ≤ −3.0 V, the device remains at the HRS, decreasing the T J and the v values.

Conclusions
In this work, the 2-D simulation of the RS processes of Si-NCs/SiO 2 ML has been analyzed. The Si-NCs were assumed as agglomerates of fixed V O which are not affected for deoxidation-oxidation process. The semiempirical model of VCM is simplified using the mathematical expression based on the Mott hopping model to calculate the resistive state considering the V O C. This approximation allows for the use of a different conduction mechanism. It was found that P-F and SCLC-F fit to the experimental HRS and LRS for both 1-BL and 3-BL devices. The simulation of SiO 2 /Si-NCs ML demonstrates that the CFs are preferentially formed near the Si-NCs enhancing the LRS (ON). The simulated voltage fits well to experimental I-V data to achieve the FORMING and SET. An IRS is observed in the FORMING process of the 3-BL device, which is explained by the preferential generation of V O in the sites that form the complete CFs, through Si-NCs. The deoxidation-oxidation model shows that the RESET process depends on the oxygen ion velocity since it is affected by Joule heating.